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1  Abstract 

This  report  summarizes  the  research  done  under  the  in-house  effort  titled  "Optimal 
Methods  for  Modulation  Classification."  The  purpose  of  this  research  is  to  revisit  the 
problem  of  blind  demodulation  and  develop  classification  algorithms  for  wide  set  of 
signals  types.  The  effort  explores  two  methodologies  in  Decision  Theory:  Likelihood  Ratio 
Test  and  inference  methods  using  Expectation  Maximization.  The  focus  will  be  mainly  in 
types  such  as  QAM,  MPSK  and  BPSK  Spread  Spectrum.  The  algorithms  must  have  medium 
complexity  and  provide  some  optimal  classification. 

In  this  report,  we  explored  the  challenges  of  constructing  LRT  rules  for  BPSK  spread 
spectrum.  Simple  assumptions  such  as  assuming  statistically  independent  waveform 
coefficients  proved  to  be  inappropriate  to  develop  such  rule.  An  assumption  of  using  a 
probability  distribution  based  on  the  Total  Square  Correlation  resulted  to  be  effective, 
although  it  turn  out  to  be  a  complicated  procedure  for  waveforms  coefficients  derived 
from  4x4  Hadamard  matrices. 

The  second  approach  used  was  based  on  Information  Divergence.  This  approach  can  be 
adapted  for  multiple  modulation  schemes.  It  was  effective  for  simple  modulation  schemes 
such  as  BPSK,  QPSK  and  QAM16,  but  it  has  convergence  problems  for  complex 
constellations.  Future  work  will  apply  the  same  technique  to  spread  spectrum  signals  and 
will  try  to  improve  the  convergence. 

Finally,  this  research  explored  MIMO  codes  design.  The  purpose  was  to  design  codes  that 
exceeded  the  performance  of  the  Full-Rate,  Full  Diversity  Golden  Codes.  The  approach 
consisted  in  searching  for  linear  2x2  MIMO  codes  that  can  allocate  6  symbols  instead  of  4. 
Although  it  was  possible  to  find  such  codes;  however,  none  of  the  codes  that  were  found 
exceeded  the  performance  of  the  Golden  Code. 
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2  Introduction 


Existing  methods  in  modulation  classification  can  be  grouped  in  three  categories:  heuristic 
methods,  feature  based  and  decision  theory  classification. 

Heuristic  methods  are  derived  from  the  observation  of  the  signal  in  time  and  frequency 
domains.  Examples  of  these  are:  the  Phase  Histogram  classifier  was  proposed  by  Lifdtke 
(Lifdtke,  1984)  and  the  well-known  Mth  Law  Classifier  for  MPSK  signals.  The  performance 
of  these  methods  is  unpredictable.  From  a  research  point  of  view,  heuristic  methods  have 
minor  importance  in  the  area  of  modulation  classification. 

Feature-based  methods  rely  on  the  heuristic  selection  of  features,  but  introduce  a  more 
scientific  approach  by  using  clustering  algorithms  for  decision  making.  These  methods  are 
more  popular  in  the  literature  and  they  make  use  of  statistical  moments.  (Dobre,  Abdi, 
Bar-Ness,  &  Su,  2007)  They  use  data  set  for  the  characterization  of  the  signals  under 
classification.  The  performance  of  these  methods  is  known  to  depend  on  the  training  data. 
Examples  of  these  methods  are  the  Statistical-Moment  Based  Classifier  presented  by 
Soliman  (Soliman  &  Sue,  1992)  and  the  Sombrero  Classifier  developed  by  Army  RDECOM. 

Finally,  the  decision  theoretic  methods  are  based  in  Bayes  Decision  Theory.  These 
methods  offer  superior  performance  because  they  are  optimized  for  minimum  probability 
of  error.  (Huang  &  Polydoros,  1995)  Methods  based  on  optimal  decision  theory  are  less 
abundant  in  the  literature  due  to  the  complexity  in  the  development  of  such  algorithms. 
They  require  good  models  and  some  appropriate  optimization  criterion  for  deriving  the 
optimal  detector.  The  algorithms  are  designed  to  guarantee  the  best  performance  in  low 
SNR.  The  few  known  methods  of  this  kind  are  signal-specific  and  most  of  them  apply 
simple  MPSK  and  QAM  models. 

The  purpose  of  this  research  is  to  extend  the  decision  theoretic  methods  to  broader  sets  of 
signals  while  trying  to  minimize  the  complexity  of  the  algorithms.  To  achieve  this  goal,  we 
will  explore  two  paths:  the  development  of  likelihood  approaches  and  the  use  of  an 
entropy  based  method  known  as  the  Minimum  Divergence  Principle  (MDP). 

The  discussion  of  this  report  is  organized  as  follows:  First,  the  MPSK  Likelihood  Ratio  Test 
will  be  revisited  and  modified  to  accommodate  spread  spectrum  signals.  As  an  alternative 
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approach,  inference  methods  such  as  EM  and  MDP  will  be  applied  to  MPSK/QAM  and  then 
extended  to  spread  spectrum.  Apart  from  decision  theoretic  methods,  the  report  explores 
the  design  of  MIMO  codes  for  optimal  performance.  MIMO  codes  are  considered  a  spatial- 
temporal  modulation  with  the  purpose  of  understanding  MIMO  codes  and  extending 
classification  algorithms  to  this  set  of  signals. 
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3  Modulation  Classification  Methods 
3.1  Existing  Likelihood  Methods 

Various  classification  methods  using  decision  theory  can  be  found  in  the  literature.  These 
methods  targets  a  specific  set  of  signal  that  include  a  combination  of  MPSK  and  QAM. 
Table  1  shows  some  of  these  developments.  (Dobre,  Abdi,  Bar-Ness,  &  Su,  2007) 

Table  1.  Existing  Decision  Theoretic  Methods 


Method 

Signal  Types 

Authors 

ALRT 

BPSK,  QPSK 

Huang,  Polydoros  (Kim  & 

Polydoros,  1988) 

ALRT 

MPSK 

Huang,  Polydoros  (Huang  & 

Polydoros,  1995) 

ALRT 

16PSK,  16  QAM 

Sapiano  (Sapiano  &  Martin,  1996) 

GLRT 

16PSK,  16  QAM 

Panagiotou  (Panagioutou, 

Anastasoupoulos,  &  Polydoros, 

2000) 

HLRT 

16PSK,  16QAM,  64QAM 

Dobre  (Dobre,  Modulation 

Classification  in  Fading  Channels 

using  Antenna  Arrays,  2004) 

EM 

BPSK-SS 

Yao,  Poor  (Yao  &  Poor,  2000) 

The  list  is  not  exhaustive,  but  helps  to  understand  some  of  the  challenges  in  the  area  of 
research.  Average  Likelihood  Ratio  Test  (ALRT)  requires  taking  the  expectation  of  all 
unknown  variables.  The  unknown  variables  are:  the  symbols,  which  are  usually  several 
hundreds  or  thousands;  the  phase  offset;  and  the  time  offset.  Integrating  over  hundreds 
symbols  is  possible  if  the  symbols  are  considered  statically  independent,  then  the 
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expectation  (i.e.,  integration)  of  products  becomes  the  product  of  the  expectation  (i.e., 
integration). 

Some  developments  have  tried  to  simplify  the  integration  by  performing  a  finite  sum. 
These  methods  are  known  as  Quasi-ALRT.  Other  methods  just  provide  a  good  estimate  of 
some  unknown  parameters.  The  methods  are  known  as  Generalized  LRT.  Other  methods 
contain  a  mixed  of  methodologies.  These  methods  are  called  Hybrid-LRT. 


The  methods  in  Table  1  may  correspond  to  one  of  the  types  listed  in  Table  2. 

Table  2.  Summary  of  Modulation  Classification  Methods 


Likelihood  Based 

Feature  Based 

Likelihood  Ratio  Test  (LRT) 

Expectation  Maximization 

Approaches  (EM/MDP) 

Statistical,  Histogram  and 

Spectral  Features 

Goal:  Classification  of  MPSK 

or  QAM 

Goal:  Classification  of  MPSK,  QAM 

and  Spread  Spectrum 

Goal:  Classification  of  MPSK, 

QAM 

Approach:  Derivation  of  a 

rule  for  classification  using 

the  signal  model 

Approach:  Inference  of  the 

likelihood  distribution  front  the 

correlator  output 

Approach:  Generation  of 

features  using  peak  power  and 

frequency  information; 

classify  signals  using 

deviations  in  frequency  or 

amplitude 

Assumptions:  Test  a  single 

test  signal, 

Cases:  Coherent  Detection, 

Non-coherent  Detection,  or 

Asynchronous  detection 

Assumptions:  Test  a  single  signal, 

Requires  time  synchronization, 

Non-coherent  Detection 

Assumptions:  Test  multiple 

non-overlapping  signals  in  a 

wide  spectrum.  May  have 

multiple  modulation  types. 
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Likelihood  Based 


Feature  Based 


Likelihood  Ratio  Test  (LRT) 


Expectation  Maximization  Statistical,  Histogram  and 

Approaches  (EM/MDP)  Spectral  Features 


Signal  Model: 


>(0  =  I  V2£  bk  eje°p{t 


-  kT ) 


Stochastic  Model: 


p(r|z;  9) 


(tiNo)1/2 


Signal  Model: 


>(n)  =  ^  Sk(n)ej 


j(2nfrn+0rni 


+  N(ri) 


r(t)  =  s(t)  +  n(t) 
n(t)~AWGA 


M— 1 

=  y _ l_ 

/-I  M(nN0 ) 


r  exp  — 


j2nl 

r  —  e  m 


Mixture  of  signals  is  possible 


Classification  Criteria: 


Classification  Criteria: 


Classification  Criteria: 


Likelihood  ratio  test: 


EpU (riff,)} 

(r)  EpU(r|H0)}^ 


Conditional  Likelihood: 


Minimum  Divergence 


Dkl  (q(r,z;  0)||  p(r,z;  0)) 


After  parameter  estimation: 


Decision  Tree  Classification 
using  features  front  histogram 


X  —  eNo^°  Relr(t)s  (t.piHilldt  =  arg  max Z)KL  (q(r,z;0)||  p(r,z;0) 


Strengths: 

Optimal  performance 


Strengths: 

Does  not  require  complex 
equations,  extended  to  Spread 
Spectrum 


Strengths: 

Easier  implementation 


Weaknesses: 


Weaknesses: 


Weaknesses: 


Complex  decision  rules  or 
non-analytical.  Limited 
implementation  to  MPSK  and 


Convergence  problems;  multiple 
local  minima 


Susceptible  to  power 
fluctuations;  Inability  to 
differentiate  narrow  band 
signals  and  tones. 


Approved  for  Public  Release;  Distribution  Unlimited. 
6 


3.2  Likelihood  Ratio  Test  (LRT) 

One  of  the  most  relevant  papers  found  in  modulation  classification  is  "Likelihood  Methods 
for  MPSK  Classification".  This  was  the  first  method  in  departing  from  feature-based 
classification  by  developing  simple  rules  for  classifying  the  set  of  MPSK,  i.e.,  BPSK,  QPSK, 
8PSK,  16PSK  and  others.  The  method  validated  existing  heuristic  Mth  power  law.  It 
showed  that  the  developed  rules  were  special  cases  of  a  more  general  decision  rule.  The 
method  served  as  a  starting  point  for  developing  classifiers  of  spread  spectrum  signals. 

3.2.1  MPSK  Model 

We  consider  a  baseband,  finite  sequence  of  symbols  bk  with  a  phase  offset  e7®0  carried  by 
non-overlapping  normalized  pulses  p(t).  The  number  of  symbols  is  K  and  the  duration  of 
each  pulse  is  T.  The  signal  s(t)  is  a  superposition  of  pulses  modulated  by  the  respective 
symbols,  as  shown  in  (1).  The  energy  of  each  symbol  is  represented  by  E. 

K-l 

s(t)  =  ^  VE  bk  e-,e°p(t  —  kT)  for  0  <  t  <  NT  ^ 

k= 0 

For  MPSK  information  symbols  the  symbols  have  amplitude  of  one  and  discrete  phases 
which  are  uniformly  distributed  around  the  unit  circle,  (2).  The  M  represents  the  number 
of  symbols.  In  the  case  of  BPSK  signal,  M  =2. 

bk=eie*,  =  k  =  0,1 . M  (2) 

3.2.2  Decision  Theoretic  Approach 

The  likelihood  function  is  just  a  density  function  conditioned  to  a  given  hypothesis.  (3)  If 
the  stochastic  model  is  a  composite  model  with  several  parameters  beside  the  random 
variables,  the  probability  density  function  is  referred  as  a  conditional  likelihood  function. 

KX\ H)  =  P(X\H )  (3) 

The  likelihood  ratio  is  the  ratio  of  two  likelihood  functions  under  two  hypotheses. 
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(4) 


A(X)  = 


KX\Ha) 

KX\Hb) 


When  we  deal  with  signals,  we  talk  about  processes.  A  process  is  a  random  event  that  can 
be  characterized  with  a  parametric  probability  density  function  (pdf).  The  random 
variable  X  is  replaced  with  all  the  possible  indexed  variables  r(t).  If  r(t)  is  a  set  of 
Gaussian  independent  variables  the  variance  is  infinite: 

tf 

A(r(t)|H)  =  lim  FT  -=e-b(t)-s(t)l2/«o  f51 

At^0  U 

t=t+At 


Redefining  a  likelihood  ratio  for  process  requires  introducing  the  notion  of  orthogonal 
basis  representation.  A  process  n(t)  will  be  decomposed  in  arbitrary  but  deterministic 
basis  functions  as  shown  in  (6).  Using  the  assumption  of  a  wide  sense  stationary 
process  (7)-(8)  and  treating  the  transmitted  signal  as  a  deterministic  process,  it  is  possible 
to  compute  the  likelihood  ratio  for  a  finite  bandwidth  process.  (See  Appendix  A) 


nit)  =  ^n^(t) 


1  =  0 


(6) 


E[n(t )}  =  0 


(7) 


N0 

E{n(t)n(u )}  —  —S(t  —  u) 


(8) 


For  an  AGWN  process,  the  likelihood  functions  are  defined  according  to  (9).  The 
hypothesis  H  is  implicitly  stated  in  the  model  chosen  for  the  transmitted  signal  s(t). 


A(r(t)|tf)  =  e 


(r(t)-s(t))2 

N0 


(9) 


Most  of  the  complexity  of  (9)  appears  when  we  consider  unknown  parameters  such  as 
symbols  bk,  phase  offset  60  (phase  synchronization),  symbol  period,  and  time  offset  (time 
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synchronization).  In  this  case,  the  conditional  likelihood  needs  to  be  averaged  over  all 
possible  symbols. 

3.2.3  Case  1:  Averaging  Over  Unknown  Symbols 

Averaging  over  symbols  is  done  under  the  assumption  of  a  known  distribution.  In  the 
worst  case,  it  is  assumed  that  the  symbols  are  equally  distributed.  In  the  case  of  MPSK,  we 
have: 


p(bk  =  eie)  =  ^  8{bk  -  e'e) 

(10) 

M  — 1 

V{bk\MPSK)  =  -  2^  8{bk  ~ 

(11) 

m= 0 

f  (r(t)-s(t))2) 

A(r(t) \H)  =  Ebobl^bK\e  N°  j 

(12) 

The  statistical  independence  of  each  symbol  allows  separating  the  expectation  of  all 
symbols  into  the  product  of  the  expectations  (13).  This  may  bring  some  computational 
benefits. 


IE 


&0&1  1 


(r(t)-s(t))2 

N0 


(r(t)-s(t))2> 

N0 


1=0 


(13) 


3.2.4  Case  2:  Averaging  Over  Unknown  Phase 

The  previous  equation  is  a  case  of  the  GLRT  when  the  phase  offset  90  is  unknown.  The 
generalized  likelihood  method  uses  the  best  estimate  of  90  rather  than  averaging  over  all 
possible  values  of  the  parameter. 

For  MPSK  case,  the  distribution  of  90  is  assumed  to  be  uniform,  i.e.,  there  is  no  preference 
in  the  values. 


pW  =  ^ 


(14) 
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(15) 


A(r(t)|H)  =  E0l 


E 


b0b1...bK 


(r(Q-s(fr))2>|l 

Ij 


3.2.5  Case  3:  Averaging  Over  Unknown  Timing 

Perhaps,  the  most  important  parameters  are  time  offset  and  symbol  period.  Recovering 
the  symbols  would  be  almost  impossible  without  any  of  these  parameters.  Yet,  many 
literature  papers  make  the  assumptions  that  unknown  signals  are  perfectly  synchronized. 
One  may  argue  that  a  preamble  must  be  present  in  any  communication  system  and  thus  it 
is  possible  to  synchronize  signals  with  high  degree  of  accuracy. 

The  time  offset  synchronization  is  modeled  by  an  error  in  the  limits  of  the  integration  in 
the  correlator  as  shown  in  (16)  and  (17).  Equation  (18)  is  the  overall  expression  for  the 
likelihood  of  a  signal  with  unknown  data,  phase  offset  and  time  offset. 


r(k  +  l)T 

rk(e)  =  I  r(t  -  eT)p(t)*dt 

JkT 


(16) 


r(k+l-e)T 

rk(e)=  r(t)p(t  +  eT)*dt 

J{k-e)T 


(17) 


(  (r(Q-s(t))2> 

=  EAEg*Ebobl'„bK  \e 


(18) 


3.2.6  Development  of  MPSK  Classifiers 

The  development  of  optimal  decision  rules  tries  to  simplify  the  equations  usually  by 
changing  the  order  of  the  expectations,  reducing  the  equations  using  Taylor 
approximation  and  using  some  other  desirable  means  of  simplification  such  as  taking  the 
log-likelihood. 

3.2. 6. 1  Averaging  Over  the  Symbols 

For  a  synchronous  MPSK  detector,  the  likelihood  function  is  given  by  (19).  A  new 
parameter  rk  is  defined  as  the  correlator  output  (20).  The  parameter  y  is  the  signal-to- 
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noise  ratio  at  the  output  of  the  correlator.  The  details  of  this  development  are  shown  in 
Appendix  B. 

K-l 

A(r(t)|H)  =  Y[Ebk{e^Re(r^bKe~Je°l-YN}  (19) 

k= 0 

2  r(k+l-e)T 

rk{e)  =  -  r(t)  p(t  -  (fc  -  e)T)  dt  (20) 

v  J(k-e)T 

Averaging  over  an  arbitrary  symbol  provides  the  standard  form  of  the  likelihood  function. 

K- 1  M/2-1 

A(r(t)|H  =  M)  =  e~NKr  cosh  ^(c)  e  m  e  -,e°)j  (21) 

fc=0  m=0 

Equation  (21)  can  be  conditioned  on  phase  and  time  offset  parameters  by  taking  the 
expectation  over  those  parameters  as  shown  in  (18). 

Huang  and  Polydoros  (Huang  &  Polydoros,  1995)  noted  in  a  binary  hypothesis  test  of 
A(r(t)|M)  and  A(r(t)|M'),  where  M  >  M' and  both  M  and  M'  are  powers  of  two,  it  was 
possible  to  reduce  (20)  into  two  expressions.  The  first  expression  is  a  common  term 
shared  by  A(r(t)|M)  and  A(r(t)|M')  as  a  function  of  M'  and  the  second  term  is  a  function 
of  M.  The  term  provided  in  (22)  will  serve  as  the  decision  rule  for  MPSK  signals. 


qM  =  ln{Mr(t)\M )}  -  Zn{A(r(t)|M)} 


The  term  (rfe(e))M  is  the  heuristic  power  law  detector  for  MPSK  signals.  Given  an  ideal 
MPSK  signal  r(t)  =  eJ2na(k)/Mi  the  power  law  detector  produces  a  constant  response 
r(t)M  =  eJ2naW  in  the  baseband  signal  and  a  fixed  tone  in  the  cosine  modulated  signal. 

The  rule  for  the  non-coherent  detection  case  is  shown  in  (22).  Later  on,  during  the 
discussion  of  the  development  of  the  spread  spectrum  detector,  we  will  show  that  the 
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detection  rule  is  identical.  The  conclusion  is:  it  is  not  possible  to  distinguish  a  spread 
spectrum  signal  from  an  MPSK  signal  based  on  simple  assumptions  of  independent 
waveform  coefficients. 

3.2. 6.2  BPSK  Spread  Spectrum  Model 

The  case  of  spread  spectrum  signals  shows  some  similarities  to  the  MPSK  classifier.  In  this 
scenario,  we  deal  with  complex  waveforms  instead  of  simple  pulses.  The  waveforms 
belong  to  a  set  of  mutually  orthogonal  or  quasi-orthogonal  waveforms  according  to  (23) 
and  (24).  Each  waveform  will  serve  as  a  separable  coding  channel.  The  spread  spectrum 
signal  is  the  superposition  of  modulated  waveforms  according  to  (25). 

L—l 

pu(t)  =  —  'YjCuiipit  -  (i  -  1)T)  0  <t<TL  (23) 

*  1  i= 0 


(24) 


s(t)  = 


\ 


U-1K-1 


L—l 


LTeJd°  X  X  buk  X  °ui  ^  “  kLT  ~  ^  ~  ^  ^ 

u= 0 k=0  i= 0 


(25) 


Under  MPSK  hypothesis  with  unknown  phase  and  time  offsets,  the  likelihood  takes  the 
following  form: 


A(r(t)|  H) 


=  EAEf 


C  j  b00b01-bU-l,K-l  I  C00c01-cU-l,L-l  I' 


\e~k^TS  Re[  r(t)-s*(C-e)}dt-j^/0'V7’s|s(t-e)lzdt ) 


(26) 


The  energy  term  of  the  received  signals  was  removed  from  the  equation  because  it  is 
constant.  The  transmit  signal  term  is  kept  because  the  amplitude  of  the  spread  spectrum 
signal  is  no  longer  constant.  It  is  easy  to  note  the  complexity  of  this  expression.  The 
expectation  must  be  taken  over  each  possible  combination  of  waveform  coefficients  and 
symbols.  In  the  case  of  binary  spread  spectrum  signals,  the  integral  would  require  2KL  +  2 
integrations. 
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We  started  a  simple  development  by  averaging  over  simplistic  assumptions.  We  begin 
assuming  that  the  waveform  coefficients  and  symbols  are  statistically  independent.  This 
would  allow  having  an  equation  similar  to  (19),  where  the  expectation  of  product  equals 
the  product  of  expectations. 

3.3  Development  of  LRT  for  BPSK-SS 

3.3.1  Simple  Averaging  Over  the  Waveform  Coefficients  and  Symbols 

A  synchronous  BPSK  spread  spectrum  detector,  the  likelihood  function  is  given  by  (27).  A 
new  parameter  rk  is  defined  as  the  correlator  output  (28).  The  parameter  y  is  the  signal- 
to-noise  ratio  at  the  output  of  the  correlator.  The  details  of  this  development  are  shown  in 
Appendix  B. 


A(r(t)\H)  —  E, 


c00c01  —CU 


(  —f 

\eN0  Jo 
-l.L-1  ( 


Re{r(t)-s*(t-e)}dt-^-fgTs\s(t-£)\2dt 


) 


(27) 


2  /-(fc+e)T 

rki(e )  =  —r =====  r(t)  xp(t  -  (k  -  e)TL  -  (i  -  1)  T)  dt 

J NnLT, 


^N0LT( 


c  J(k-l+e)T 


(28) 


Assumption  #1:  The  energy  term  of  the  transmitted  signal  will  be  treated  as  a  global 
constant.  The  distribution  of  the  overall  symbol  amplitude  (29),  which  is  the  sum  of  the 
symbols  times  the  waveform  coefficients,  is  a  binomial  distribution. 


U-1K-1  L—l 


;[m]  =  ^  ^  buk  ^  cui  S(m  -  kL  -  (i  -  1)  ) 


(29) 


u= 0 k= 0  i= 0 


The  term  is  E, 


c00c01 


■cu-i,l-i  | 


-2- dVTs|s(t-e)|2dt)  .  r  .  £  ,  ■  rrni-i 

e  No  0  is  a  function  of  the  variance  of  s  [mj  which 


depends  on  the  number  of  users  U  and  length  of  the  spreading  code  L.  So  the  only  term  of 
interest  is  the  correlation  term. 

Assumption  #2:  The  waveforms  coefficients  are  statistically  independent.  The  likelihood  is 
simplified  to  (30). 
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(30) 


U-l  L- 1 


Mr(t)\H)  ±Y\Y\Ec*t  {eW°J° 


u= 0  i= 0 


■^7/0NTsRe{r(t)-s*(t-6)}dt| 


It  was  found  that  the  assumption  is  not  useful.  Averaging  over  statistically  independent 
symbols  will  also  yield  a  product  of  terms  as  in  (31). 


u- 1  K-l 


fU-l  L- 1 


A(r(t)|H)  =  f]  f]  fff)  ECul  ( 


u'= 0 fc  =  0 


^u=0  i= 0 


(31) 


A  lengthy  development  reveals  that  using  a  second  assumption  results  is  equivalent  to  the 
BPSK  likelihood.  The  conclusion  is:  the  assumption  of  statistically  independent  coefficients 
does  not  allow  differentiating  between  a  BPSK  signal  and  a  BPSK  spread  spectrum  signal. 

3.3.2  Weighted  Average  over  the  Waveform  Coefficients  and  Symbols 

The  proposed  approach  uses  a  statistically  dependent  distribution  for  averaging  over  the 
waveform  coefficients.  The  approach  requires  us  to  design  a  waveform  density  based  on 
some  metric  that  characterizes  spreading  sequences.  We  consider  the  metric:  Total  Square 
Correlation  (TSC)  defined  as  the  Frobenius  norm  of  the  correlation  matrix. 

Let  C  —  {cui}  be  the  matrix  of  waveform  coefficients  row  u  and  column  i.  Also  consider  a 
row  vector  representation  cu  —  [cu0,  cul, cul-i\-  The  correlation  matrix  is  defined  as 
(32).  (Cotae,  2003) 


u- 1  u- 1 


TSC(C)  =  ^  ^  \cuc%\2  >U,  for  L  >  U 


(32) 


u= o  u'= o 


A  modification  of  this  metric  was  used  for  our  density  as  shown  in  (33). 


u-i  u- 1 


u= 0  u  —  0 


u±ur 


(33) 
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For  a  binary  spread  spectrum,  the  coefficients  take  two  real  values:  cui  —  ±1,  so  any 
power  also  stays  in  the  same  range:  c” j  —  +1.  The  density  for  statistically  dependent 
coefficients  is  (34). 


e-SD(C ) 


(34) 


U-l  U- 1  L-l  L—l 


D(C )  =  (U2-U)L  +  4  1111  Cik.CilCjk.Cjl 

i=0  j= 0  k= 0  1=0 


J*i 


k=tl 


(35) 


The  expression  of  the  likelihood  (36)  is  complicated  and  further  analysis  will  be  provided 
in  a  future  report.  In  the  present  discussion,  we  are  going  to  develop  the  likelihood 
function  for  the  simplest  case  U  —  L  —  2. 


A(r(t)|H) 


=  Re{rui(.e)b*ukc*uie  ie°}-Ylu=o^i= 


kcui^ukcui  0{cuj 


coo  coi  CU-1,L-1 


(36) 


3.3.3  Simple  Case:  2x2  Hadamard  Coefficients 

The  likelihood  for  detecting  2x2  Hadamard  coded  spread  spectrum  was  computed  in  (37) 
assuming  statistical  dependency  of  the  waveform  coefficients.  The  results  were 
encouraging  as  we  can  see  from  Figures  Figure  1-Figure  3.  A  derivation  of  the  likelihood 
was  done  in  Mathematica.  (See  Appendix  C) 
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\{r(t)\H:  {  Cij  =  ±1,U  =  2,L  =  2}) 

=  4.8.(i  +  e2g)  (e~2^  +  e"25cosh  (2V7  Re(rkoj) 
+  e~25cosh  [2^Jy  Re(rkl)^j 
+  e2S~4r  (cosh  (2  Jy  i?e(rfc0))  cosh  (2  Jy  Re(rkSjfj 


(37) 


s  =  V 27 


Log-Likelhood  BPSK  @  5  dB 

600  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - r 


150  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

0  200  400  600  800  1000  1200  1400  1600  1800  2000 

test  label 


Figure  1.  Likelihood  of  BPSK 

The  likelihood  of  the  spread  spectrum  signal  can  be  decomposed  in  two  terms:  the 
individual  BPSK  terms  for  rk0  and  rkl>  and  the  product  of  both  terms.  The  LRT  test  gives 
acceptable  results  for  SNR  ranges  between  0  and  5  dB.  For  higher  ranges,  the  test  has 
problems  due  to  overflow  (1000  chips  used)  and  the  test  breaks  for  levels  below  -5dB. 
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Figure  2.  Likelihood  of  BPSK  Spread  Spectrum 


Figure  1  and  Figure  2  show  the  likelihood  for  BPSK  and  BPSK  Spread  Spectrum.  Both 
likelihoods  can  be  approximated  using  Gaussian  distributions.  The  empirical  receiver 
operating  characteristic  (ROC)  curve  is  shown  in  Figure  3. 


ROC  Curve 

H0:BPSKvs.H1:BPSK-SS 


Figure  3.  ROC  Curve  for  LRT  test 
HO:  BPSK,  HI:  BPSK-Spread  Spectrum 
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3.4  Inference  Methods  for  Classification 

In  order  to  overcome  the  complexity  of  our  blind  demodulation  problem,  this  paper 
proposes  to  modify  a  stochastic  approach  for  signal  demodulation  by  incorporating  the 
concept  of  information  divergence  or  "distance"  between  probability  density  functions. 
One  of  the  densities  provides  a  parametric  model;  the  second  density  describes  the  ideal 
modulation  type.  (See  equations  (38,  (56)  and  (58))  The  models  can  describe  noise  or  any 
arbitrary  modulation  type.  An  optimization  algorithm  (54)  will  adjust  the  parameters  of 
the  distribution  to  fit  various  models  by  minimizing  the  divergence.  The  decision  is  made 
by  selecting  the  smallest  "distance"  instead  of  the  likelihood  ratio. 

The  Expectation  Maximization  approach  proposed  by  Yao  (Yao  &  Poor,  2000)  provides  an 
iterative  method  for  estimating  the  parameters  in  equation  (38).  The  scheme  uses  BPSK- 
SS  model  with  a  known  number  of  users.  The  signal  model  is: 

y  —  Vfi7  ■  S  -b  +  n  (38) 

Where  n  is  the  noise,  b  E  {+1}*-  is  the  transmitted  symbol  vector  of  length  K.  The  matrix 
S  is  a  L  x  K  matrix  with  K  users  and  L  the  signature  size.  The  S  matrix  contains  the 
normalized  spreading  sequences  S  —  [s0,  sv ...%]  E  is  the  energy  of  the  signal.  The 
output  y  is  the  received  vector.  The  process  can  be  described  with  a  correlation  matrix  R : 

E  {yyH}  =  R  =  UHAU  (39) 

Eigenvalue  decomposition  provides  a  means  to  separate  the  signal  space  from  the  noise 
space.  The  highest  eigenvalues  correspond  to  the  energy  of  the  signals  and  the  associated 
eigenvectors  correspond  to  the  signature  vectors  that  form  matrix  S. 

A  =  As  +  An  =  diag(A0,  Xlf ... ,  AK,  AK+1, ... ,  AL) 

As  =  diag(A0,A1,  ...,AK)  signal  space  (40) 

An  =  diag(  Ak+1,  ... ,  Al)  noise  space 

The  associated  eigenvectors  Us  to  the  signal  space  can  be  used  to  estimate  the  matrix  S  by 
means  of  (41). 
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U"y  —  x  —  Ve  ■  S  -b  +  o  ■  n 


(41) 


The  conditional  likelihood  function  becomes  (42): 

2k 

p(x|VIs)  =  exp(__Lp_VssSt||2)  (42) 

k. — 1 

If  bk  is  a  binary  vector,  then  there  are  2K  possible  combinations  of  symbols.  The  EM 
process  adds  a  new  variable  to  the  conditional  likelihood.  This  parameter  is  known  as  the 
missing  data  in.  The  missing  data  is  a  label  that  identifies  the  population  or  cluster  where 
the  known  variables  belong  to.  In  this  case,  the  cluster  is  any  binary  pattern  formed  from 
numbers  0  to  2K  —  1. 

N  i 

p(x,in|VfS)  =  ^ - Iexp(-^p-  VIShjQpOJ  (43) 

n=i{2nos)2  s 

The  EM  method  maximizes  the  expectation  of  the  log-likelihood  given  the  data.  This  is: 

S  =  arg  max  Einly  {log  (p(x,in  | VfTS))}  (44) 

A  detailed  explanation  of  this  algorithm  explained  in  (Yao  &  Poor,  2000).  By  using  an  EM 
approach,  the  process  of  maximizing  the  expectation  eventually  that  converges  to  the 
waveform  coefficients  in  matrix  S  without  dealing  with  the  complexity  of  the  ALRT. 
According  to  the  paper,  the  estimation  of  the  waveform  coefficient  offers  a  slightly  worse 
performance  compared  to  a  minimum  square  error  estimator  during  the  demodulation. 
The  worse  performance  is  expected  because  MMSE  assumes  known  waveform  parameters 
while  in  the  EM  approach  the  estimation  is  totally  blind. 

The  development  of  the  EM  estimator  allows  for  applying  constraints,  but  these  are  not 
explored  in  the  referenced  paper.  For  instance,  other  possible  constraint  is  to  minimize 
the  entropy  of  the  likelihood  function  (45)  for  all  possible  signals.  This  constraint  would 
force  that  the  results  reduce  the  randomness. 
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(45) 


N 

^p(x|in)log(p(x|in))  =  0 

n= 1 

3.4.1  Information  Divergence  Modulation  Classifier 

Under  the  FY11-12  AFOSR  Mini-Grant,  the  concept  of  information  divergence  was  studied 
as  a  symbol  estimator  of  unknown  signals.  The  concept  was  used  as  a  blind  estimator  as 
reference  (Ding  &  Kay,  2011)  suggested. 

Instead  of  using  a  ratio  of  likelihood  functions,  the  proposed  approach  uses  the  Kullback- 
Leibler  (KL)  divergence.  KL  divergence  is  a  measure  of  "distance"  or  dissimilarity  between 
two  distributions.  Given  two  distributions  u(x,  y)  and  v(x,  y),  the  Kullback-Leibler 
divergence  is  defined  by  (46): 


DKl(u\\v)  =  -  j  J  u(x,y)ln  dxd7  C46) 

If  two  distributions  are  the  same,  then  it  is  easy  to  observe  that: 

Dkl(v\\v)  =  0  (47) 

The  divergence  between  two  distributions  is  always  positive.  (Cover  &  Thomas,  2006) 

Dkl(u\\v)>0  (48) 

The  divergence  is  related  to  the  likelihood  ratio.  Consider  two  distributions  v(x,  y)  and 
u(x,  y),  then  the  divergence  between  v  and  u  is: 

DKl(u  lb)  =  Eu(r,Hl){/op(  A(r))}  +  Eu(r,Hl)  j  log  (j^y)}  C49) 


where  A(r)  is  the  likelihood  ratio: 


/l(r)  = 


u(r  jtfi) 

v(r\H0) 


(50) 
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and  Eu(r>H  )  is  the  expectation  operator.  Using  divergence  is  an  extension  of  the  concepts 
found  in  Decision  Theory.  The  divergence  includes  an  average  of  the  logarithm  of  the 
likelihood  ratio. 

KL  divergence  has  been  used  in  the  speaker  recognition  area.  In  this  problem,  a  feature 
space  vector  x  contains  information  that  is  used  to  discriminate  between  speakers. 
Suppose  that  one  has  a  stochastic  model  p0(x)  such  as  a  multivariate  Gaussian 
distribution  J\f(p0,£0)  using  the  features  of  the  population  as  random  variables.  Suppose 
that  a  second  Gaussian  model  describes  a  target  speaker.  A  way  for  designing  a  decision 
rule  Pi(x)  =  N(jl1,'Z1)  between  the  two  classes  consists  of  using  a  subset  of  xthat 
maximize  the  divergence  DKL(p1\\p2)  between  the  two  classes.  The  distance  between  the 
data  and  the  model  is  computed  using  only  the  selected  features  and  the  class  that  has  the 
"shortest  distance"  wins.  The  concept  is  similar  to  the  Mahalanobis  distance  d l,  — 
(x  —  pj)TEtrl(x  —  /fj)  used  to  discriminate  between  classes.  In  fact,  the  KL  divergence  is 
equivalent  to  the  Mahalanobis  distance  when  p0  =  ^  for  unimodal  Gaussian  distributions 
(Campbell,  1997). 

The  proposed  modulation  classifier  uses  similar  concepts  of  minimum  distance  for 
deciding  between  classes.  The  feature  space  is  the  complex  plane  that  contains  the  values 
of  the  estimated  symbols,  each  one  with  an  in-phase  and  quadrature  component.  The 
stochastic  model  is  assumed  to  be  Gaussian.  The  true  symbols  form  a  constellation  which 
can  be  rotated  due  to  phase  incoherency  between  the  transmitter  and  the  receiver. 

We  define  our  communication  channel  using  the  diagram  of  Figure  4.  (Cover  &  Thomas, 
2006)  The  symbols  can  be  drawn  from  any  arbitrary  distribution.  It  can  be  circular  or 
rectangular  QAM  or  it  can  have  arbitrary  constellations. 


Figure  3.  Communication  Channel 
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The  probability  of  a  symbol  p(z)  is  assumed  to  be  an  unknown  variable  for  computational 
purposes.  This  prior  knowledge  must  match  the  implicit  probability  of  symbols  in  the 
template,  as  it  will  be  explained  later. 


Table  3.  Definitions  for  EM  Approach 


Concept: 

Mathematical  Expression: 

Symbol 

z  G  Szwith  Sz  being  a  finite 

discrete  set  of  complex  numbers 

Symbol  Probability 

P(z) 

Observation 

r  E  Sr  with  Sr  being  a 

finite  set  of  complex  numbers 

Noise  Power  Density 

2  _  No 
°  2 

Likelihood  of  the  observation 

given  a  symbol  (channel  noise) 

p(r|z)  ~  M(z,  <r2) 

Probability  of  a  symbol 

given  the  observation 

p(r|z)p(z) 

P(z  r)  = 

P(r) 

Matching  density  function 

q(r) 

Sampled  Based 

Kulback-Leibler  Divergence 

Dkl(u\\v)-  ££u(r,z)h.(  ) 

TESr  ZESZ  V  ' 

MAP  (Maximum  Posterior  Rule)  for 

choosing  between  classes  Cl  and  C2 

Ct 

p(z|r,C1)^p(z|r,C2) 

Co 
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The  channel  introduces  noise  of  each  symbol.  The  noise  model  is  expressed  as  a  likelihood 
function.  The  noise  is  assumed  to  be  AWGN  but  this  model  can  be  modified  to  other  types 
of  noise  if  needed. 

The  correlator  in  the  receiver  produces  the  observation  r.  For  each  observed  data  point 
there  is  a  probability  of  a  symbol  given  the  observation  p(z|r).  The  maximum  posterior 
probability  criterion  can  be  used  as  a  decision  rule  to  compute  the  detected  symbol  z. 

Bayes’  theorem  states  that: 


p(z|r)p(r)  =  p(r|z)p(z)  (51) 

Note  that  the  divergence  between  p(z|r)p(r)  and  p(r|z)p(z)  is  zero: 

Ari(p(z|r)p(r)||  p(r|z)p(z))  =  0  (52) 

Suppose  now  that  we  do  not  want  to  use  p(z).  Instead,  we  would  like  to  use  a  template 
distribution  q(r)  .  This  template  satisfies  the  inequality: 

Ari(p(z|r)<jf(r)||  p(r|z)p(z))  >  0  (53) 

The  equality  holds  only  and  only  if  q(r )  =  p(r;  0)  for  a  given  parameter  vector  0.  If  the 
equality  does  not  hold,  there  is  at  least  a  minimum  divergence  that  can  be  achieved 
between  the  two  distributions  according  to  Kullback’s  Minimum  Discrimination  of 
Information  Principle  (MDI).  (Gray,  2001)  The  parameter  vector  can  include  noise,  phase 
and  signal  energy  parameters.  In  the  simulations,  only  phase  incoherency  was  considered. 

The  optimal  parameters  can  be  determined  using  the  proposed  optimization  criterion: 

§°Pt  _  arg  max  Dkl  (q(ji  z;  p)  ||  p(r>  z;  fl))  (54) 

e 

where  q(r,  z;  0)  =  p(z|r;  0)q(r)  and  p(r,z;  0)  =  p(r|z;  0)p(z;  0).  The  algorithm  can  be 
regarded  as  a  stochastic  matched  filter. 
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Once  the  optimal  parameters  have  been  determined  for  classes  with  templates  (?i(r)  and 
q2(j),  the  decision  is  made  either  using  the  conventional  maximum  posterior  probability 
or  using  the  KL  divergence. 

3.4.2  Simulations 

MPSK  models  and  QAM  models  were  considered  for  the  estimation  of  the  parameters. 

3. 4.2.1  MPSK  Case: 

If  symbols  z  belong  to  the  MPSK  class,  then: 

(  y  2  71771 

z  G  je  m  ,m  =  0, 1,  ...,M  —  lj  (55) 

The  likelihood  model  that  describes  the  channel  is  the  simple  Gaussian  distribution  shown 
in  (56): 


,  .  1  /  |r  —  ze>B  I  \ 

p(rM)  =  o^CJtp( - —)  (56) 

The  model  describes  a  noisy  output  with  rotated  symbols  z  e;e.  The  parameter  vector 
consists  of  two  quantities,  the  phase  mismatch  and  the  probability  of  the  symbol: 

9  =  [9,  p(z)]  (57) 

The  template  of  an  arbitrary  MPSK  hypothesis  is: 

q(r\HM)  = 

Note  that  the  template  incorporates  the  true  prior  distribution  p(z)  =  1/M  .  This  is  not  a 
necessary  requirement,  but  provided  a  convenient  quick  implementation. 

Figure  4  to  Figure  8  shows  the  convergence  of  the  algorithm  (49)  for  MPSK  signals.  The 
observed  data  r  is  shown  as  green  dots.  The  red  dots  track  the  convergence  of  the  data  to 


M—l 

I 

1=0 


M(_nN0y/2 


exp 


j2nl 

r  —  e  m 


N  n 


(58) 
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the  model,  i.e.,  ze]  -*  e  m  .  The  blue  dots  represent  the  true  symbols  and  the  yellow  dots 
(barely  visible)  represent  the  final  convergence.  The  black  lines  represent  the  decision 
boundaries  for  MPSK  types.  (See  Figure  5  and  Figure  6) 


d 


Minimum  Divergence  2PSK  Model  end  2PSK  Dalai 
180  Iterations.  SNR  =  3  dB 
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Figure  4.  Convergence  of  the  MDI  Estimator/Classifier 
using  BPSK  Data 


In  Figure  5,  we  achieved  convergence  in  4  iterations.  BPSK  is  the  simplest  form  of  MPSK 
types.  Quick  convergence  is  expected  even  in  low  SNR  levels. 


Approved  for  Public  Release;  Distribution  Unlimited. 
25 


15 


Minim  um  Diverg  enc  e  Q  P  S  K  M  o  d  el  an  d  Data 
180  Iterations,  SNR=  5  dB 


10 


5 


■s 


o 


-5 


-10 


•15  - 1 - 1 - 1 - l 

-15  -10  -5  0  5 

Inpha3e 


15 


Figure  5.  Convergence  of  the  MDI  Estimator/Classifier 
using  QPSK  Data 


In  Figure  5,  convergence  of  QPSK  data  to  QPSK  model  is  achieved  in  approximately  50 
iterations.  The  computational  time  is  still  insignificant  due  to  the  implementation  of  the 
algorithm  which  is  based  on  simple  linear  algebra  and  element  by  element  operations. 

In  Figure  6,  convergence  of  QPSK  data  to  BPSK  model  is  achieved  in  approximately  50 
iterations.  The  minimum  distance  between  a  QPSK  and  a  BPSK  densities  is  achieved  when 
we  distribute  the  data  equally  along  the  decision  line.  The  resulting  divergence  of  this 
case  is  greater  than  the  divergence.  In  such  case,  the  classifier  would  decide  for  QPSK 
instead  of  BPSK. 
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Figure  6.  Convergence  of  the  MDI  Estimator/Classifier 
using  QPSK  Data  and  BPSK  model 


The  reader  might  wonder  if  an  algorithm  such  as  Gaussian  Mixture  Models  (GMM) 
clustering  algorithm  can  perform  a  similar  task.  The  answer  is  negative.  Differently  from 
the  GMM  algorithm,  the  MDI  classifier  keeps  fixed  proportions  between  the  symbols  (or 
clusters)  and  prohibits  them  from  moving  independently. 

Figure  7  shows  a  case  of  a  noisy  16PSK  signal.  The  convergence  of  cases  like  this  one  may 
take  several  hundred  of  iterations  but  running  in  fractions  of  seconds.  It  is  very  likely  that 
optimal  classifiers  fail  in  a  case  like  this  one.  The  only  solution  to  a  problem  like  this  is  to 
observe  more  data. 
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Minim  um  Diverg  ence  for  1 6MPSK  lype  1 80  iterations  1 5  dB 


Figure  7.  Convergence  of  the  MDI  Estimator/Classifier 
using  16PSK  Data 


3.4.3  QAM  Case 

The  proposed  classifier  can  be  flexible  enough  to  incorporate  other  constellations  such  as 
QAM  models.  (Figure  8)  The  only  requirement  is  to  provide  a  template  as  an  input.  Other 
constellations  (templates)  are  admissible  by  just  providing  a  different  template 
model  q(r|H). 

The  output  of  the  algorithm  is  the  optimal  parameter  vector  and  the  divergence.  With  the 
parameter  vector,  the  likelihoods  and  priors  for  each  observation  can  be  reconstructed. 
This  is  good  feature  because  these  variables  not  only  allow  deciding  the  modulation  type, 
buy  they  provide  sufficient  information  for  recovering  the  symbols  and  the  bitstream. 

In  the  case  of  complex  constellations  (QAM-32,  QAM-64),  the  convergence  is  less  reliable 
due  to  the  multiple  local  minima  that  exist.  The  method  is  more  likely  to  work  on  simple 
constellations  like  BPSK,  QPSK  signals  and  BPSK-SS  waveforms. 
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Figure  8.  Convergence  of  the  MDI  Estimator/Classifier 
using  16QAM  Data 


3.4.3. 1  Divergence  for  Waveform  Estimation 

The  benefit  of  using  information  divergence  is  that  the  method  preserves  the  signal 
constellation  proportions  and  only  allows  rotation  or  scaling.  This  is  a  feature  the  EM 
method  does  not  have. 
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4  MIMO  Code  Design 

The  purpose  of  this  research  study  is  to  investigate  the  MIMO  Space-Time  Block  Codes 
(STBC)  and  provide  if  possible  an  improved  design.  In  this  report,  the  available  codes  will 
be  presented.  The  best-known  code  will  be  identified  and  study  for  the  purpose  of 
developing  better  codes. 

4.1  Introduction 

MIMO  is  a  multiple  input  multiple  output  communication  systems  that  transmit 
information  from  a  multiple  antenna  transmitter  to  a  multiple  antenna  receiver  in  the 
presence  of  fading.  One  of  the  benefits  that  make  MIMO  system  attractive  is  the  increasing 
channel  capacity  due  to  the  antenna  diversity.  MIMO  systems  are  currently  used  in  mobile 
networks  such  as  cell  phones  and  routers,  where  voice  and  data  requires  higher 
transmissions  rates  and  better  quality  of  service. 

The  design  of  In  MIMO  codes  is  one  important  topic  of  MIMO  systems.  The  codes  found  in 
literature  can  be  classified  in  the  groups  as  follows: 

■  Layered  Spaced  Time  Coding 

■  Space  Time  Block  Coding 

■  Trellis  Space  Time  Coding. 

Each  code  family  differs  in  the  degree  of  complexity  versus  accuracy.  Layered  Space  Time 
Coding  is  an  example  of  simplistic  codes  with  a  poor  performance  compared  to  trellis 
codes  with  the  highest  complexity  and  the  best  accuracy.  Block  coding  provides  a  trade-off 
between  performance  and  complexity. 

This  report  is  dedicated  exclusively  to  the  study  of  the  STBC.  The  design  of  STBC  does  not 
provide  error  correction  as  the  Trellis  Codes,  only  mitigation  of  multipath.  The  best- 
known  STBC  is  known  as  the  Golden  Code.  It  was  developed  around  2003  and  since  this 
date  no  other  code  with  superior  performance  has  been  found. 

Desirable  characteristics  that  make  a  good  code  can  be  described  in  terms  of  the  following 
properties: 

■  Best  coding  gain 
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■  Minimum  average  power 

■  Full  rate 

■  Lowest  complexity 

After  a  literature  search,  the  research  was  oriented  to  study  and  propose  a  new  scheme 
that  exceeds  the  performance  of  the  Golden  Code.  For  such  purpose,  we  will  present  in  the 
discussion  the  idea  behind  the  development  of  the  Golden  Code  before  proposing  schemes 
and  testing  the  possible  advantages  or  disadvantages  of  such  codes.  Then,  we  will  proceed 
to  investigate  schemes  that  can  provide  the  properties  of  good  codes.  Finally  we  will  make 
an  assessment  of  the  designs  and  get  to  conclusions  whether  there  is  a  better  STBC  than 
the  Golden  Code. 

4.2  Space-Time  Model 

The  problem  of  finding  new  codes  will  be  constrained  to  the  case  of  frequency  non- 
selective  channels,  also  known  as  the  narrowband  case.  We  also  consider  a  two-transmit 
two-receive  antenna  system  in  the  entire  discussion  unless  it  is  otherwise  specified.  In  this 
scenario,  the  frequency  response  of  the  channel  coefficients  is  flat  in  the  bandwidth  of 
interest.  The  channel  matrix  in  (59-(60)  characterizes  the  channel  fading: 


/in  ^12 

-h-2i  h22. 


(59) 


where  /ii;  is  the  gain  of  the  path  between  transmitter  antenna  j  and  receiver  antenna  /.  The 
relationship  between  the  output  and  the  input  signals  under  fading  and  additive  white 
Gaussian  noise  is:  Y  —  H  ■  C  +  N  or, 

phi  yi2i  _  pin  ^12 1  _  rcn  ci2i  phi  ni2i 

bhi  y 22 J  [h21  h22\  lc2 1  c22\  |_n2i  H22J 

where  yi;-  is  the  received  signal  i  at  time  ci;-  is  the  transmitted  symbol  from  receive 
antenna  /  at  time  and  is  the  noise  seen  by  receiver  antenna  i  at  time  j. 
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4.2.1  Criteria  for  Best  Performance 

4.2. 1.1  Coding  Gain 

In  MIMO  systems,  the  role  of  the  code  is  to  transmit  information  in  such  way  that  it 
provides  immunity  to  fading  in  AWGN.  A  metric  of  performance  is  the  minimization  of  the 
pairwise  probability  of  error  (PEP).  The  PEP  given  a  set  of  transmitted  codes  Ck  and  a 
detected  code  Ck  is  defined  as: 


Pe 


k= 0 
k*j 


(61) 


The  analysis  (Paulraj,  2003)  shows  that  the  PEP  a  function  of  the  rank  r,  the  number  of 
receive  antennas  Mr,  the  product  of  the  eigenvalues  Ak of  the  difference  Cj  —  Ck.  An  upper 
bound  for  the  PEP  is  given  by: 


Based  on  PEP  analysis  the  code  gain  is  given  by  (63): 


J 


(62) 


(63) 


For  code  design,  we  consider  a  fixed  physical  system  were  the  number  of  transmit  antenna 
Mr  equals  the  number  of  receive  antennas  Mr.  We  also  considered  the  case  when  the 
number  of  time  periods  equals  the  number  of  transmit  antennas  such  as  the  one  shown  in 
(60). 

Optimizing  PEP  requires  maximizing  the  code  gain.  This  is  achieved  in  two  ways:  the 
matrix  difference  Cj  —  Ck  is  full  rank.  This  is  the  rank  criterion.  The  second  is  the 
determinant  of  the  matrix  is  as  large  as  possible.  The  combination  of  the  two  criteria 
yields: 
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J  =  min 


(64) 


4.2. 1.2  Minimum  Average  Power 

The  design  the  best  space-time  code  requires  finding  a  set  of  codes  Ck  such  that  the  the 
set  maximizes  the  coding  gain  and  minimizes  the  average  power.  We  look  at  the  following 
Frobenius  norm  or  a  matrix  trace  for  optimality. 


max  \det{Cj  -  Cfc)| 
min  rr(Cfc) 


(65) 


For  comparison  purposes  between  codes  sets  A  and  B  with  the  same  code  rate,  we 
normalize  the  power  of  A  and  B  and  find  the  highest  code  gain. 

4.2. 1.3  Code  Rate 

Code  rate  is  defined  as  the  number  of  symbols  that  can  be  transmitted  in  a  STBC  code.  For 
a  2x2  matrix,  the  maximum  number  of  symbols  is  4.  If  the  code  matrix  can  transmit  the 
maximum  number  of  symbols,  then  the  code  is  known  as  full  rate  FR.  If  in  addition  the 
code  has  full  rank,  the  matrix  is  referred  as  full-rate,  full  diversity  code  (FRFD).  The  best 
performance  STBC  codes  are  FRFD. 

4.2. 1.4  Code  Complexity 

A  tradeoff  between  performance  and  complexity  is  preferable  when  the  demodulation  cost 
of  the  receiver  exceeds  the  benefit  of  overcoming  the  fading  effects.  Such  cost  can  be 
measured  in  terms  of  system  power  consumption,  size  and/or  processing  time.  Under  our 
discussion,  code  complexity  was  not  a  consideration  for  the  design  of  the  best  STBC. 

4.3  Existing  MIMO  Codes 

The  first  task  in  the  search  for  the  best  performance  codes  is  to  consider  the  existing  code 
design  through  the  literature.  As  we  mentioned  early,  MIMO  codes  found  in  literature  can 
be  classified  in  the  groups. 
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4.3.1  Layered  Space  Time  Coding 


One  of  the  simplest  coding  schemes  is  called  Horizontal  Encoding  (HE).  In  this  scheme  the 
transmitted  symbols  are  demultiplexed  into  separate  streams  or  layers  through  each 
transmit  antenna  in  single  or  multiple  time  slots.  Vertical  Encoding  (VE)  is  another 
variant  where  the  same  symbol  is  spread  through  all  antennas  as  opposed  of  being  spread 
through  multiple  time  slots.  Diagonal  Encoding  (DE)  transmits  the  symbols  using  a 
different  antenna  at  a  different  time  slot.  (Sellathurai  &  Haykin,  2009)  Codes  like  VE  or  DE 
are  in  general  easy  to  implement,  but  not  accurate  due  to  error  propagation. 

4.3.2  Space  Time  Block  Coding 

Similar  to  Layered  Space  Time  Codes,  these  blocks  repeats  symbols  and  its  complex 
conjugates  in  a  way  that  the  code  achieves  full  diversity. 

4.3.2.1  Alamouti  Codes 

Alamouti  proposed  a  2x2  coding  scheme  that  meets  full  diversity  and  simple  decoding. 
(Blahut,  2003)  The  code  rate  is  one  symbol  per  time  slot.  The  design  proposed  design  was 
extended  higher  dimensionalities.  Its  code  form  can  be  expressed  as: 


Alamouti  codes  are  linear  codes  that  are  orthogonal  according  to: 

C  ■  CH  =  (\Xl\2  +  \x2\2)  I  (67) 

This  property  facilitates  the  simplification  of  the  maximum  likelihood  demodulator: 

C  =  arg  min||V  —  H  •  C\\\  (68) 


C  —  arg  min  —2 


Tr{Y  ■  H 


H-C}-^{\x1\2  +  \x2\^\\H-Hh\\22 


(69) 


Further  simplification  allows  expressing  the  demodulator  by  expressing  the  demodulator 
in  separable  functions  ofxi  andx2. 


Approved  for  Public  Release;  Distribution  Unlimited. 
34 


Mathematically  speaking,  the  Alamouti  code  is  a  Hamiltonian  Quaternion  isomorphism.  A 
Hamiltonian  Quaternion  is  an  algebraic  structure  known  as  ring  that  obeys  associative  and 
distributive  laws  under  addition  and  multiplication  while  preserving  the  same  algebraic 
structure  of  the  operands.  This  discovery  was  found  latter,  with  the  development  of  FRFD 
codes. 

4.3.2.2  Quasi-Orthogonal  Codes 

Several  variants  of  Alamouti  code  for  higher  dimensions  have  been  proposed  (Sharna  & 
Papadias,  2003).  These  code  matrices  are  not  orthogonal  but  the  demodulator  can  be 
expressed  in  separable  functions.  For  higher  dimensionality  the  code  rate  of  an  Alamouti 
and  quasi-orthogonal  matrices  decrease,  so  the  transmission  becomes  less  efficient. 

4.3. 2. 3  Algebraic  Codes 

These  codes  are  developed  based  on  algebraic  number  theory.  The  goal  of  these  codes  is  to 
ensure  full  diversity  and  full  rate. 

4.3. 2.4  Golden  Code 

Golden  Codes  were  presented  by  (Belfiore  &  Viterbo,  2004)  as  cyclic  division  algebra,  but 
other  isomorphic  forms  of  the  code  were  already  known.  The  Golden  Code  can  transmit  4 
symbols  in  two  time  slots  and  ensure  full-diversity.  So  far,  it  is  claimed  that  the  Golden 
Code  has  the  best  performance  of  all  known  codes. 

x0  X1  I 

ya(x0)  ffOi)]  j-70-j 

x0  —  a  +  9b  xx  —  c  +  9d 

4.3.2. 5  Perfect  Space-Time  Codes 

Perfect  codes  are  an  extension  of  the  2x2  Golden  Codes  for  higher  dimensions. 

x0  xx 
ya(x0)  a(x±) 
y2cr(x0)  yoixi) 

-Y Mx0)  y2o-(*i) 


*2  *3 

ct(x2)  er(x3) 
cr(x2)  cr(x3) 
ya(x2)  a(x3). 


(71) 


Approved  for  Public  Release;  Distribution  Unlimited. 
35 


It  has  been  said  that  cyclic  algebra  codes  provides  full-rate  linear  codes  while  cyclic 
division  algebras  provide  full  diversity  codes. 

4.3. 2. 6  Trellis  Codes 

Trellis  codes  offer  the  best  performance  at  the  expense  of  a  complex  implementation. 
They  are  also  referred  as  convolutional  codes.  Due  to  the  complexity  of  the  design,  these 
codes  are  excluded  from  this  report. 

4.4  Group  Theory  and  Golden  Codes 

The  concept  of  an  algebraic  structure  describes  a  numeric  set  such  as  numbers,  vectors, 
polynomials,  matrices  of  an  invariant  form  that  defines  one  or  more  operations  such  like 
addition,  subtraction,  multiplication  and  division. 

Groups  are  algebraic  structure  under  addition  and  subtraction.  The  operands  belong  to 
the  same  domain  and  the  result  stays  in  the  domain  and  preserves  the  same  form  of  its 
operands.  The  operands  obey  associative,  distributive  and  identity  laws  under 
addition/subtraction  or  multiplication.  Examples  of  groups  are  2x2  matrices  defined  in 
the  real  or  complex  domain. 


CL  -y 
Cl 


d-2. 


a3  b3 
-c3  d3 


CLi,  bi,  Ci,  di  G  TL 


(72) 


Rings  are  algebraic  structures  under  addition,  subtraction  and  multiplication.  As  in 
groups,  the  operands  belong  to  the  same  domain  and  the  result  stays  in  the  domain  and 
preserves  the  same  form  of  its  operands.  It  rings  obey  to  associative,  distributive  and 
identity  laws  under  addition,  subtraction  and  multiplication.  Examples  of  rings  are 
quotient  rings  of  polynomials,  modular  arithmetic  and  triangular  matrices. 


b1 

cl2  b2 

a3  b3 

.  0  d1 

.  0  d2. 

.  0  <^3- 

di,  bi,  Ci,  dt  E  R 


(73) 
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Ideals  are  subsets  of  rings  with  special  properties.  The  product  of  any  element  of  a  ring 
and  ideals  results  in  an  ideal.  Examples  of  ideals  are  the  set  of  even  numbers.  Any  real 
number  multiplied  by  an  even  number  becomes  an  even  number. 

/  =  TL  =  (2  •  a|a  G  Z) 

a  ■  b  —  c  (74) 

if  a  E  I,  b  E  TL,  then  c  ETL 

Fields  are  algebraic  structures  under  addition,  subtraction,  multiplication  and  division.  As 
in  all  the  other  algebraic  structures,  the  result  of  an  operation  stays  in  the  domain  of  the 
operands. 

4.4.1  Isomorphism  and  Homomorphism 

Algebraic  structures  can  be  mapped  into  other  algebraic  structures.  We  talk  about 
isomorphism  when  there  is  a  one-to-one  transformation  between  structures. 
Homomorphism  is  similar  to  isomorphism  except  that  several  element  of  an  algebraic 
structure  can  be  mapped  into  a  single  element  of  the  second  structure. 

4.4.2  Why  Group  Theory? 

In  STBC  design,  we  look  for  matrices  with  full  rank  and  full  rate.  Full  rank  implies  the 
existence  of  a  multiplicative  inverse.  In  addition  to  this,  we  require  that  the  difference 
between  any  two  coding  matrices  is  non-singular.  These  two  features  can  be  found  in 
rings. 

4.4.2.1  Example  ofAlamouti  Code 

Alamouti  codes  are  an  isomorphism  of  Hamiltonian  quaternions.  The  algebraic  structure 
defines  a  basis  of  one  real  and  three  complex  quantities  {  1,  I,  J,  K  }  that  obeys  to  the 
following  laws: 


I  ■  J  =  -J  ■  I ,J  ■  K  =  I ,  K  ■  I  =  J  and  I  ■  J  ■  K  =  - 1  (75) 

Any  quaternion  can  be  express  as  linear  combination  of  these  four  group  elements: 

q  —  %  +  a2  ■  I  +  a3  ■  /  +  a4  ■  K  (76) 
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The  sum  of  quaternions  is  also  a  quaternion: 


Qi  =  a4  +  a2  •  /  +  a3  •  J  +  a4  ■  K 

q2  =  b1  +  b2-I  +  b3-J  +  b4-K  (77) 

<?1  +  <?2  =  al  +  ^1  +  (a2  +  ^2)  ‘  J  +  (a3  +  ^3)  ‘  /  +  (a4  +  ^4)  '  K 

The  product  of  a  quaternion  is  a  quaternion: 


q±  ■  q2  =  OA  -  a2h2  -  a3b3  -  a4b4 )  +  (a^  +  a2hx  +  a3b4  -  a4b3 )  ■  / 


+  (a1h3  —  a2b4  +  a3b!  +  a4b2 )  ■/ 


+  ( a4b4  +  a2b3  —  a3b2  +  a4b4 )  ■  K 


(78) 


A  quaternion  conjugate  and  inverse  is  defined  as: 

q4  —  a4  —  a2  ■  I  —  a3  ■  J  —  a4  ■  K 


such  that 


q4  ■  q4  —  af  +  af  +  a3  +  a4. 


(79) 


(80) 


An  important  isomorphism  of  the  quaternion  is  the  2x2  matrices. 

qt  =  (a4  +  a2  ■  /)  +  (a3  +  a4  ■  7)  ■  /  (81) 


Substituting  base  /  by  the  imaginary  number  i  one  can  map  quaternions  into  matrices: 


<7i  - 


a4  +  a2i 
—a3  +  a4i 


a3  +  a4i 
a4  +  a2i 


<?2  G 2 


'  b4  +  b2i 

b3  +  b4i' 

yi 

y2] 

—b3  +  b4i 

hi  +  b2i . 

L-yi 

y2J 

(82) 


There  is  a  mapping  between  addition  and  multiplication  operations  that  can  be  verified. 
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^1  +  Cl  +  C2 


(83) 


qi-q2^Gi-  G2 

A  study  of  Group  Theory  revealed  that  Alamouti  matrices  are  an  isomorphic  group  that 
obeys  the  same  rules  of  complex  quaternions. 

4A.2.2  Extension  of  Algebraic  Structures 

Extending  an  algebraic  structure  means  that  there  is  a  broader  structure  such  as  a  group, 
ring  or  field  that  includes  the  elements  of  the  structure.  Consider  the  case  of  the  set  of  real 
integers  Z.  The  set  forms  a  ring.  This  field  can  be  extended  to  a  more  general  class  that 
includes  real  and  complex  integers  denoted  as  Z[i\  also  known  as  Gaussian  Numbers. 

TL  TL[i\  (84) 

An  element  of  the  extended  structure  can  be  represented  using  the  basis  {  1,  i  }.  Any 
number  in  the  ring  can  be  expressed  as  a  pair  of  two  real  coefficients: 


x  —  a  +  b  ■  i 


x  —  a  +  b  ■  i  —  (a,  b),  a,bER 


(85) 


A  rational  number  is  defined  as  a  number  that  can  be  expressed  in  terms  of  two  integers:  a 
numerator  and  a  denominator. 


Q  —  { a/b  |  a,b  G  TL }  (86) 

The  set  of  rational  numbers  forms  a  field.  The  field  can  be  extended  to  the  domain  of 
complex  integer  numbers. 


Q  =  {a/b  | a,b  E  Z[i]}  (87) 

The  extended  field  holds  four  real  integers.  It  can  be  extended  further  to  include  the  set  of 
algebraic  numbers.  Algebraic  numbers  are  defined  as  any  possible  root  of  any  arbitrary 
rational  polynomial. 
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Q  =  {x  | P(x)  =  0 ,  P(x)  rational  polynomial } 


(88) 


A  subset  of  the  algebraic  numbers  is  known  as  the  set  of  quadratic  numbers. 

Q[i]  =  (a  +  Vd  ■  b,  a,  b  £  Q[i]}  (89) 

A  quadratic  number  can  "store"  eight  integers.  Thus  by  extending  the  field  one  can  "store" 
more  integer  coefficients  in  a  single  number  structure.  In  the  case  of  the  Golden  Code,  a 
quaternion  is  build  using  two  quadratic  numbers.  Thus,  the  Golden  Code  can  "store"  16 
integer  numbers  in  a  2x2  matrix  structure.  This  allows  getting  full  rate  and  full  diversity 
code. 

The  matrix  in  (70)  forms  a  division  algebra,  i.e.,  a  non-commutative  division  field.  The 
matrix  is  isomorphic  with  the  cyclic  algebra  defined  as: 


A  —  L±  +  e  ■  L2 


a  ■  e  —  e  ■  er(a) 


(90) 


The  value  ymust  be  a  non-relative  norm  and  the  function  o(x)  defines  commutation  rules 
for  multiplication. 

4.4.3  Algebraic  Tools  for  Code  Design 

In  the  literature,  two  main  theorems  have  been  use  in  the  design  of  MIMO  codes:  the 
Lindemann  Theorem  and  algebraic  number  theorem  found  in  (Perlis,  1952),  (Paulraj, 
2003)  and  (Ma  &  Giannakis,  2003). 

4.4.3. 1  Lindemann  Theorem 

Given  any  distinct  algebraic  numbers  ai,...am  and  non-zero  algebraic  numbers  ai,...,  am, 
then  the  following  inequality  holds: 

m 

0  (91) 

i= 1 
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4A.3.2  Decomposition  of  the  Extension  Field 

For  a  field  generator  6  in  field  F,  any  element  /?  can  be  uniquely  expressed  as  a  polynomial 
in  ^with  a,  coefficients  in  F 


m 

YJai-ei  =  p  (92) 

i= 1 

This  theorems  implies  that  element  /?=0,  then  a,-0,  but  ffe 0  implies  that  there  is  at  least 
one  coefficient  a,  that  is  non-zero. 

These  theorems  look  similar,  but  they  are  mathematically  different.  The  first  theorem  is 
based  in  the  distinction  between  algebraic  and  transcendental  numbers.  A  transcendental 
number  cannot  be  expressed  as  the  root  of  rational  polynomials.  The  second  theorem  is 
strictly  based  in  algebraic  numbers. 

4.4.3. 3  Available  Tools  for  Code  Construction 

In  the  search  for  mathematical  tools,  one  should  be  aware  of  the  existence  of  isomorphism 
of  Golden  Code.  Those  isomorphisms  may  include  similarity  transformations,  rotations, 
scaling,  linear  mappings  and  permutations.  These  forms  do  not  constitute  a  new  code 
improvement. 

4.4.3.4  Matric  Polynomials 

Matric  polynomials  of  degree  m  may  take  the  form: 

m 

^ai-Ai=A  (93) 

i= 1 

As  long  the  matrix  A  is  full  rank  and  the  symbols  a,  are  not  equal  to  the  coefficients  of  the 
characteristic  equation,  the  determinant  of  A  will  not  be  zero.  A  careful  selection  of  a 
constellation  lattice  (or  equivalently  coefficients  a,)  may  result  in  matrices  with  non¬ 
vanishing  determinants.  This  construction  was  not  explored  in  this  study. 
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4.4.3. 5  Companion  Matrices 


Companion  Matrices  are  related  to  matric  polynomials.  Code  constructions  with  non¬ 
vanishing  determinants  may  be  of  the  form: 


4 


(94) 


The  design  of  C  must  be  such  that  the  difference  between  any  code  has  non-vanishing 
determinants.  This  construction  was  not  explored  in  this  study. 

4.4.3. 6  Variations  of  the  Golden  Code 

Another  explored  form  was  matrices  with  determinants  similar  to  the  Golden  Code.  For 
the  Golden  Code,  the  determinant  is  given  by: 


det(C)  =  /(x1,x2)  -  i  ■  f  (x3,  x4) 
f(a,b )  =  (xl  xxx2  x2) 
i  ■  f(c,  d )  =  (xf  -  x3x4  -  x|) 


(95) 


We  note  that  the  determinant  has  two  quadratic  multinomials  of  the  same  form  and  one  of 
them  is  rotated  90  degrees  (multipled  by  the  imaginary  constant)  in  the  complex  plane. 
Consider  the  domain  of  a  and  b  to  be  the  set  of  Gaussian  integers  in  some  rectangular 
configuration  as  shown  in  Figure  9.  If  we  examine  the  set  of  points  generated  by 
evaluating  the  multinomial  /(a,  b )  using  all  admissible  values  of  a  and  b,  we  will  note  that 
the  co-domain  are  Gaussian  integers  and  there  are  some  empty  spaces  corresponding  to 
Gaussian  integers  that  never  occurs.  (The  codomain  is  shown  as  blue  points  in  Figure  10) 
We  also  notice  that  the  multinomial  codomain  of  i  ■  f(c,  d )  falls  exactly  in  these  empty 
places  as  shown  in  Figure  10.  The  minimum  value  of  the  determinant  of  (82)  is 
determined  by  the  minimum  distance  between  the  blue  and  red  points.  This  observation 
will  be  used  in  the  search  of  new  codes. 
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GoldenCode 
Minimum  Determinant 

20  - 1 - 1 - 1 - 1 - 


15 


>0  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

-20  -15  -10  -5  0  5  10  15  20 

Real 


Figure  9.  Complex  Plane  Showing  Co-Domain  of f(a,b) 
in  Blue  and  i*f(c,d)  in  Red 


4.5  Search  for  New  Codes 

An  interesting  question  in  the  code  design  is:  is  it  possible  to  find  codes  similar  to  Golden 
Codes  that  accommodate  more  than  four  symbols?  The  answer  turned  to  be  affirmative  in 
this  research.  Unfortunately,  the  codes  that  were  found  have  an  average  power  that  is 
several  times  higher  than  the  Golden  Code. 

The  code  structure  of  the  new  proposed  codes  is  similar  to  (95),  but  it  can  accommodate 
six  symbols  as  shown  in: 


det(C)  = 


+  a2s2  +  a3s3 

KPlS4  +  /?2S5  +  /?3S6) 


CC1S4  +  a2s5  +  a3s6 

P1S1  +  p2S2  +  P  3S3 


det(C)  =  /(Si,  s2,  s3)  +  i  ■  /(s4,  ss,  s6) 

/(si,s2,s3)  =  ( a iSj.  +  a2s2  +  a^ip^  +  p2s2  +  p3s3 ) 

=  %s^  +  a2sf  +  a3s3  +  a12S!S2  +  +  a23s2s3 


(96) 


It  was  possible  to  find  such  multinomial  by  doing  a  random  search. 
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(97) 


/(slfs2,s3)  =  %2  +  s|  -  s|  +  (17  +  13;)S!S2  +  (13  -  13 j)s±s3 
+  (-13  +  17;)s2s3 

|Re(sfc)|  <  3  and  |/m(s/c)|  <  3 

A  plot  of  the  co-domain  of  (97)  is  shown  in  Figure  10,  where  none  of  the  rotated  points  in 
red  overlap  with  the  blue  points.  The  problem  is  that  the  multinomial  does  not  factor  as  a 
product  of  linear  multinomials.  Factoring  these  terms  results  in  a  system  of  equations  that 
is  overdetermined.  No  solution  has  been  found  for  this  particular  case.  Without  this 
factorization,  the  multinomial  is  useless  as  a  code. 


Figure  10.  Complex  Plane  Showing  Co-Domain  of f(sl,s2,s3J 
in  Blue  and  i*  f(s4,s5,s6)  in  Red 


The  next  step  would  be  to  find  a  multinomial  form  that  results  in  an  overdetermined 
system,  but  has  at  least  one  solution.  The  form  of  this  code  can  be  expressed  as: 


a  i  = 


«3  = 


Ui2  it  'sj U]  2  4nq u2 
2  a2 

a 23  ±  Va|3  ~  4a2% 
2  a9 


a2 


a2 


(98) 
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Equation  (98)  is  a  solution  obtained  by  excluding  some  possible  equations  of  the 
overdetermined  system.  Some  codes  have  been  found  by  using  conjectures  about  the  form 
of  the  solution  rather  than  using  a  searching  algorithm.  These  codes  are  shown  in  Table  4. 


*  After  normalizing  code  gain. 

The  codes  found  appear  to  be  related  to  the  golden  ratio.  The  term  ^  keep  reappearing  in 
the  codes  found.  If  there  is  such  code  that  exceeds  both  performance  and  rate  of  the 
Golden  Code,  then  such  code  must  belong  to  some  sort  of  unknown  ring  in  the  complex 
domain.  The  ring,  which  acts  as  a  place  holder  for  numbers,  would  be  capable  of  storing 
more  integer  coefficients  such  as  (99).  One  may  consider  the  ring  Ok  of  the  set  K  defined 
as: 


Approved  for  Public  Release;  Distribution  Unlimited. 
45 


(99) 


K  = 


al/m  a2/m 


f  CL 


(m- 


i  )M] 


The  goal  would  be  to  find  a  scheme  that  stores  more  integers  while  keeping  the 
determinant  and  power  comparable  to  those  ones  of  the  Golden  Code.  This  type  of 
approach  will  be  explored  as  part  of  future  work. 

4.5. 1.1  Lindemann  Theorem 

In  this  study,  we  explored  the  possibility  of  finding  codes  of  the  form: 


a  ■  ea  +  b  ■  c  ■  eY  +  d  ■  es 
i(c  •  eY  +  d  ■  es )  a  ■  ea  +  b  ■  e13 


(100) 


by  conducting  an  exhaustive  search.  The  determinant  of  the  matrix  is  a  linear 
combination  of  transcendental  numbers  and  the  Lindemann  theorem  applies. 

A  common  problem  with  all  these  codes  is  that  the  determinant  approaches  to  zero  when 
the  constellation  size  grows.  Although,  Lindemann  theorem  and  field  extension 
decomposition  deal  with  different  set  of  numbers:  transcendental  and  algebraic,  it  is 
always  possible  to  find  a  transcendental  number  sufficiently  close  to  an  algebraic  number 
(Niven,  1953).  A  number  theory  theorem  states  that  any  irrational  number  can  be 
approximated  with  infinitely  many  rational  numbers  h/k  such  that: 


h  1 
k  V5  k2 


(101) 


Finding  a  good  approximation  of  an  irrational  number  is  what  causes  the  determinant  to 
go  to  zero.  If  we  consider  only  the  diagonal  elements  of  the  proposed  structure  (100), 
then: 


det(C)  =  a2e2a  -b2e2f3  =  0 


(102) 


Solving  for  zero  yields: 


a2  _  e2^ 
b2  e2a 


(103) 
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but  under  the  assumption  that  a  and  b  are  integers,  it  follows  that: 


£  —  e2  (/?-«) 


(104) 


Within  the  approximations  of  an  irrational  number,  one  must  find  the  worst 
approximation  of  rational  numbers.  For  a  case  of  two  symbols  a  and  b,  it  is  easy  to  find 
values  that  maximizes  the  determinant.  But  for  the  case  of  four  symbols,  it  the  task  is 
extremely  difficult.  A  Matlab  simulation  shown  in  Figure  11  shows  this  fact. 


Minim  um  D  et  errninant 


lm(2*beta-2*alpha) 


Figure  11.  Minimum  Determinant  of  Code 
based  on  Lindemann  Theorem  (102) 

In  Figure  11,  the  imaginary  and  real  part  of  the  exponent  has  been  plotted.  The  real  part 
goes  from  -0.3  to  0.88  using  500  steps  (shown  in  the  figure)  and  the  imaginary  part  goes 
from  0  to  271  using  400  steps.  The  color  plot  shows  the  magnitude  of  the  determinant. 
There  is  possible  to  find  codes  (represented  in  yellow  and  red)  using  Lindemann  Theorem 
such  that  the  determinant  is  greater  or  equal  than  one.  Finding  full  rate  codes  (using 
coefficients  a,  b,  c  and  d]  is  far  more  difficult  task. 
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5  Conclusions 


Two  methods  for  classification  of  modulation  types  were  developed:  LRT  for  spread 
spectrum  and  EM  based  for  QAM  and  MPSK. 

The  first  method,  which  is  a  classical  LRT,  provided  the  desired  results  when  using  a 
modified  version  of  the  TSC  metric  for  averaging  waveform  coefficients  using  a  2x2 
Hadamard  matrix.  Mathematica  provided  over  1000  exponential  terms  for  this  simple 
case.  The  LRT  test  for  waveforms  using  a  4x4  Hadamard  matrix  is  extremely  complicated. 
As  a  future  development,  we  would  attempt  to  find  either  a  recursive  formula  based  on 
2x2  Hadamard  codes  or  deducing  some  empirical  approximation  of  the  LRT. 

The  second  classification  method  provided  a  numerical  approximation  of  the  symbols  on 
conventional  modulation  types  such  as  MPSK  and  QAM.  The  convergence  of  the  algorithm 
seems  to  be  an  issue  for  complex  QAM  constellations.  As  a  future  development,  we  will 
attempt  to  improve  convergence  using  Newton’s  method  instead  of  Gradient  Descend. 
The  method  will  also  be  extended  to  BPSK-SS  following  a  Yao’s  proposed  EM  algorithm. 
(Yao  &  Poor,  2000) 

As  part  of  the  effort,  we  also  investigated  the  design  of  codes  that  could  outperform  the 
current  best-known  code:  the  Golden  Code.  It  was  possible  to  find  finite  codes  with  higher 
rates,  but  also  with  a  higher  average  power.  Future  work  will  include  the  development  of 
algorithms  that  conduct  a  constraint  search  of  these  new  codes.  The  goal  will  be  looking 
for  codes  with  higher  spectral  rates  and  similar  average  power  to  the  Golden  Code. 
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Appendix  A 

The  development  of  the  likelihood  function  for  a  process  starts  with  an  orthogonal  series 
representation. 


S(t)=  ^bjCPiCt)  (!) 

i=0 

The  signal  n(t)  is  an  Additive  Gaussian  Noise  (AWGN)  process  with  mean  zero  and 
variance  N0/2. 


n(t )  = 


1=0 


E{Ci }  =  0 

E{acj}  = 


(2) 


A  theorem  from  decision  theory  tells  that  the  sum  of  zero  mean  Gaussian  processes  is 
Gaussian,  so  the  coefficients  ct  produced  by  the  orthogonal  basis  decomposition  are 
Gaussian  random  variables. 


r(t)  =  s(t)  +  n(t)  - 


i= 0 


(3) 


The  variables  {aj  have  mean  {hj  and  variance  ^j.  The  variables  {aj  are  also 
independent. 


E{di}  =  0 

E  {|(a£  -fy)|2}  =  y  Sij 


(4) 
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Processes  s(t),  n(t)  and  r(t)  can  be  approximated  by  selecting  a  finite  set  of  coefficient.  A 
PDF  can  be  constructed  if  the  coefficients  are  finite  N. 

A  likelihood  ratio  test  can  be  constructed  as  follows: 


A (r(t))  = 


pN(r(t)\HA) 
Pn  (r(t)|HB) 


-  e 


Li= o  N0 


(5) 


The  expression  allow  us  to  arbitrary  choose  any  number  of  coefficients.  If  we  extend  the 
number  to  infinity,  then  we  have: 


A(r(t))  =  lim 

iv  — ^  oo 


pN(r(t)\HB) 


_y N  (.r(t)-s(t\HA))2-(r(t)-s(t\HB))2 

A (r(t))  -  e  2ii=0 


(6) 


The  likelihood  function  is  defined  as  a  density  function: 

_yN  ( r(t)-s(t\HA ))2 

A(r(t))  =  e  2,1=0  wo  (7) 
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Appendix  B 

BPSK  Development 

The  hypothesis  is  implicit  in  the  selected  model.  In  the  case  MPSK  signals  we  have: 


s(t)H  —  M)  — 


K-l 

1 

k= 0 


VI  e^ke^e°  p{t  —  kT )  for  0  <  t  <  NT 


(8) 


For  unknown  types,  the  phase  of  the  symbol  9k  belongs  to  a  set  of  M  elements: 


ek  = 


(9) 


The  pulse  is  defined  as  an  ideal  square  pulse  with  period  T.  The  likelihood  can  be 
simplified  to: 


K-l 

A(r(t))  =  Y\  Ebk  {e^Reibk^idcrk^-y] 

k=0 


(10) 


where  y  is  the  SNR  at  the  output  of  the  decorrelator  and  rki(e )  is  the  decorrelator  output.  After 
some  math,  the  likelihood  expression  becomes: 

A(r(t))  =  {cosh{yJyRe{e~ie<=rkie)}))  (11) 


BPSK  Spread  Spectrum  Development 

This  development  assumes  uniformly  distributed  coefficients. 

r2n 

pbk  ieTyRe{bke-^cy,t- £  Cj*  rfci(e)}-yj  _  I  eylYRe{bke~iec  c.«  rfci(e)}-y  p(bk)dbk 

1  ’  Jo 

PiVk)  =  \  ^  5  {bk  ~  eJn^)  ^ 

m= 0 

For  q  e  M 
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(13) 


K- 1  L 

A(r(0)  =  ]~ [  ]~ [  e_r  cosh(A/7/?e{e_-,6,crfei(e)}) 

fe=0  i= 1 

The  correlation  term  rfej(e)  can  be  substituted  byrfe(e).  This  means  that  for  this  simple  case 
the  likelihood  of  BPSK  and  BPSK  spread  spectrum  are  the  same.  The  ratio  test  will  not  be 
able  to  differentiate  between  BPSK  and  spread  spectrum  under  the  assumption  of 
statistically  independent  symbols. 


Approved  for  Public  Release;  Distribution  Unlimited. 
54 


Appendix  C 

BPSK  Spread  Spectrum  Development  with  Weights 

The  equations  were  derived  using  Mathematica.  This  is  the  likelihood  when  using  a  2x2 
Hadamard  matrix. 


In [337]:=  (*  CDMA  *) 

Clear[b,  c,  r] 

M  =  2; 

H  =  Array  [Subscript  [e,  ttl  -  1,  1x2  -  1]  /  Sqrt [M]  &,  {M,  M}]; 
B  =  Array  [Subscript  [b,  xi-1,  K]  &,  {M,  1}]; 

R  =  Array  [Subscript  [r,  K,  nl  -  1]  &,  {M,  1}]; 

X  =  Rt.H.B 
Y  =  (H.B) 


Out[342}= 


* 


C0.0  *K.O  c1.0*K.l 


•  b, 


iCo.i^ic.o  Ci.xric.ij 

l  'JT  nJ~2~  ) 


}} 


Out[343)=  || 


^S.kCo.O  rb0,KCx,0  ^l.KCl 


V  2 


VT  1  VT  +  VT  ^ 


ln[344]:=  (*  WEIGHTS  *) 

TOT  =  0; 

For  [  i  =  0,  i  <  M,  i  ♦+,  For  [j  -  0 ,  j  <  M,  j  +  +  , 

For [  k  =  0,  k  <  M,  k  +  +  ,  For [1  =  0,  1  <  M,  1++, 

If  [  i  #  j  &&  k  #  1, 

TOT  +=  Subscript  [c,  i,  k]  *  Subscript  [c,  i,  1]  *  Subscript  [c,  j,  k]  *  Subscript  [c,  j,  1]  ; 

] 

n 

n 

FML  -  Expand  [  (  (M  A  2  -M)  M  +  TOT)  /  4  ] 


Out[346}=  1  +  co.o  Ca.i  Ci.o  Ci.i 


ln[347]:=  (*  ENERGY  *) 

ET  =  (YT.Y /  2)  ; 

ET  =  Expand [ET] ; 

ET  =  ET  /.  v_^2  »  1; 
ET 


Out[3SO]=  ||l 


1 

^.K^.K^D.O  Co.i 

2 


1 

2 


^0.  K  i^l.K  C1.0 


Ci.i 


}} 


ln[351]:=  (*  CONDITIONAL  LIKELIHOOD  BPSK-SS  *) 

CL  =  Exp  [Sqrt  [2  y]  •  (X  -  FML)  -  y  ET]  [[1,  1]]; 
Z  =  Join  [Sqrt  [M]  »Flatten[H],  Flatten[B]]; 
For  [  i  =  1,  i  £  Length  [Z] ,  i  ++, 
rl  =  {  Z  [  [i]  ]  -  z  ); 

CL  =  CL  /.  rl; 

CL  =  1  /  2  *  Sum  [CL,  (z,  l-l.  1)}]; 

] 

Expand [CL] 


Out[3S4)=  1  +  1  c-rtVTr,,  +  1  e-»*jV7rlt  +  1  6-Y-jV7rt:  + 

4  8  8  8 

^  -Y-2  v  Y  xK  ■  ^  “2  Y-Y  2  Vy~  !  -2-  V1' 2  ri.j-v's  rK,  i  ^  -2Y-V2  VT  f-2*V2  rx.o-vT  rK,» 

8  16  16 

-2Y*VT\‘Y  I-2-V2  rz.e^vT  r*.i|  -2  y-VT  \  y  |  -2-V  2  rs.s-VTrx.il 

16  16 
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BPSK  Spread  Spectrum  Development  with  Weights 

The  equations  were  derived  using  Mathematica.  This  is  the  likelihood  when  using  a4x4 
Hadamard  matrix. 


In(319]:' 


OuH324}- 


OutJ325}- 


lnJ326]:' 


Out(328J- 


(★  CDMA  *) 

Clear[b,  c,  r] 

M  =  4; 

H  =  Array  [Subscript  [c,  si  -  1,  s2  -  1]  /  Sqrt  [M]  &,  {M,  M}  ]  ; 
B  =  Array  [Subscript  [b,  s-l,  K]  &,  {M,  1}J; 

R  =  Array [Subscript [r,  K,  Si  -  1]  &,  {M,  1}]; 

X  =  rt.h.b 
Y  =  (H.B) 


f  f  ( 1  1  1  1 

<  bo.r  -  cc.o  *k.o  +  -  =i.e  rE1  +  -  c2.5  rE2  +  c3,c  rE,3 
I  2  2  2  2 

|  1  1  1  1 

bi,E  -  co.i  rj;,c  +  -  ci.t  rx.i  +  C2.1  rx.a  +  —  C3.1  rz.i 

I  2  2  2  2 

|  1  1  1  1 

b2,x  =0,2  rx,o  +  c3,2  rE>3  +  c2,»  rE,2 +•  c3,2  *x,3 

I  2  2  2  2 

I  1  1  1  1 

b3.x  ca.j  *x.o  +  Cl. 3  rE,3  +  c2,3  rE,2  +  c3,3  rE,3 

2  2  2  2 


}} 


•I  /  ^ 

b3.x  co,3  |  •  1  -  bo,x  =i.o  • 
2 

bi.x=2.i+-  b2,K  cJi2  *  -  bj-E  c3-3  J,  {-b0tEc3,o 


M  bo.x=e.o  +  bl-E  ca,3  *  b2,Ece,2 
1 1 2  2  2 

r  1  1 

{  b0.K  =2,0  +  - 


bi,rci,i  ♦  b2.x=i.2  +  b3E  c1>3 

2  2  2 

111 
bi,x=3,i  +  b2,Ec3,2  +  b3-Ec3,3 


(*  WEIGHTS  *) 

TOT  =  Of 

ror[  i  =  0,  i  <  M,  i++,  For [  j  =  0,  j  <  m,  j  +  +  , 

For[  Jc  =  0,  k  <  M,  k  +  +  ,  For[l  =  0,  1  <  M,  1+  +  , 

If [  i  *  j  it  k  #  1, 

TOT  +=  Subscript[c,  i,  k]  *  Subscript [c,  i,  1]  *  Subscript [c,  j,  k]  *  Subscript [c,  j,  1]  ; 

] 

]) 

]) 

FML  =  Expand  [  (  (M  '  2  -  M)  M  +  TOT)  /  4  ] 


12  ♦ 

c 

e,o  Co,i  c 

1.0  Cl.l  +  c 

0.3  C 

0,2  C 

1.0  Cl 

2  *  Co.i  Co, 2  c 

•1,1  c 

2  ♦  Co,0 

Co.i 

C 

1,0  c 

3  ♦ 

C0r 

t 

ce.s 

ci.i 

Cl.l 

♦ 

ca,2 

Co. 3 

ci.2 

=1.3 

+ 

c0,3 

=0.1 

c2,e 

c2,i 

+ 

Cl,3 

clr 

1 

c2 

0 

c2.i 

+ 

C; 

0 

ce,2 

c2.0 

c2,2 

♦ 

ci,o 

C1.2 

c2.3 

c2.2 

+ 

=0.1 

=0.2 

C2,l 

c2,2 

+ 

=1,1 

C1 

2 

c2 

1 

c2.2 

+ 

C2 

o 

ce.s 

c2.0 

c2, 3 

♦ 

ci,o 

Cl. 3 

c2.0 

=2.3 

+ 

C3,l 

=0.3 

c2 , 1 

c2,3 

+ 

=i.i 

=1 

3 

1 

C2.S 

+ 

C2 

2 

co.a 

c2.2 

C2.3 

♦ 

=1.2 

=1.3 

c2.2 

=2.3 

+ 

c3,3 

c3,i 

ca,o 

c3, 1 

+ 

=1.0 

1 

0 

c3.1 

+ 

-2 . 3 

=2.1 

=3.0  =3.1 

+ 

C5 

0 

CC.2 

c3.3 

c3.2 

♦ 

ci,0 

=1.2 

c3.0 

=3.2 

+ 

c2,3 

=2,2 

c3,0 

=3,2 

+ 

=0.1 

c3 

2 

=3 

1 

c3.2 

+ 

cl.i 

ci,2 

=3.1  =3.2 

+ 

C2 

1 

c2.2 

c3.1 

=  3.2 

+ 

c3,0 

=0.3 

=3.o 

=3.3 

+ 

=1.0 

=1.3 

ca,e 

c3,3 

+  ^2,0 

=2 

3 

=3 

2 

=3.3 

+ 

CC.l 

=  0.3 

=3.1  =3.3 

+ 

Cl 

1 

Cl.  3 

C3.1 

ca.a 

♦ 

Ci,  1 

C2.3 

C3.1 

C3.3 

+ 

C3, 2 

Co,  3 

C3.2 

C3,3 

+  Cl,  3 

Cl 

3 

C3 

2  C3,S 

+ 

C2.2 

C2,3 

C3,2  C3i3 
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lf*32S]:-  (*  ENERGY  *) 

ET  =  (yt.  Y /  2)  } 

ET  =  Expand [ET]; 

ET  =  ET/.  v_*2  1; 

ET 


Out(332)-  j|2*-  b5-K  bj.,,;  Ca,0  C0rl 

1 

—  b3  .X  c0,l  c0,3  + 

4 

1 


—  bt.x  bj.j;  Ci.i  ci, 3  + 


1 

—  be, k  bi.j;  cj,o  cj,  1  + 


1 

bl.it  b3,jr  Cj.i  Cj,3  + 


1 

—  biiK  bj,};  C3.1  C3,3  + 

4 


1 

4  4  4 

111 

—  bj.K  b3jE  cs,2  c0-3  +  —  be  E  bii}:  Ci,e  Ci,i  +  —  bc,};  bjjr  Ci,0  c2-2  * 
111 

-  b«,K  bs,r.  ci,o  Ci.j  +  —  bi,  r.  b3 ,3:  ci.i  ci,s  +  —  bt.xbj.x  ci,i  Ci.s  + 
111 

—  bs.K  b j,r.  cj,«  cj.a  +■  —  bi,r  bj.r.  cj,i  cj.a  +  —  bc.K  bs,r  cj.o  cj.s  + 
111 

-  bJtS  ba  ,K  Cj.j  Cj,3  +  —  b0iK  bi.j;  Cs,e  C3.1  +  —  bc,E  bj.i;  C3 ,  n  C33  + 


111 

-  b0.sbj,r.  C3,o  Cj,3  +  -  biiKbs.,;C3,i  Cj,3  +  -  b3,Kb3,KC3,3  C3 , 3 

4  4  4 


ln[353]:-  (★  CONDITIONAL  LIKELIHOOD  BPSK-SS  *) 

CL  =  Exp [Sqrt  [2  y]  *  (X  -  FML)  -  y  ST)  [  [1/  1)  ]  ; 
Z  =  Join [Sqrt [M]  *  Flatten [H],  Flatten[B]]; 
For [  i  =  1/  i  3  Length [Z]  ,  i  ++, 
rl  =  {  Z[  [i]  ]  -»  z  }S 
CL  =  CL  /.  rl; 

CL  =  1  /  2  *  Sum  [CL,  {z,  (-1,  1}}]; 

] 

Expand [CL] 


Out(336}- 


3  c-<>  f*  -fT 

4331 


It  + 


:  .  =  -/t 


••it 


•j’s  V» 


1,  .  1, 


2C4£ 

T 


3  «*•  V  a  -f7 

1024 


«1221» 


-•V*  V  * 


k  ;  * 


-•  r.'JT  JT 


■ti-3  i*  :.j n  -.art  »x  3 


The  likelihood  function  has  1228  exponential  terms.  The  terms  are  combinations  of  the 
vectors  components  of  the  spread  spectrum  signal. 
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List  of  Acronyms 


ALRT 

BPSK 

BPSK-SS 

Dkl 

GLRT 

HLRT 

LRT 

MIMO 

MPSK 

No 

Pdf 

PEP 

QAM 

ROC 

STBC 

TSC 

0 

R 


Average  LRT 

Binary  Shift  Keying 

BPSK  Spread  Spectrum  or  CDMA 

Kullback-Leibler  Information  Divergence 

Generalized  LRT 

Hybrid  LRT 

Likelihood  Ratio  Test 

Multiple  Input  Multiple  Output  System 

M-ary  Shift  Keying 

Normalized  Noise  Power 

Probability  Density  Function 

Pairwise  Error  Probability 

Quadrature  Amplitude  Modulation 

Receiver  Operating  Characteristic 

Space-Time  Block  Codes 

Total  Squared  Correlation 

Golden  Ratio 

Real  Numbers 


Approved  for  Public  Release;  Distribution  Unlimited. 
58 


Rational  Numbers 


1 

Complex 

m 

Gaussian 

Gaussian 

Numbers 

Integers 

Distribution  with  mean  /r  and  variance  o 2 
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